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ABSTRACT 

We make use of a semi- analytic cosmological model that includes simple prescriptions 
for dust attenuation and emission to make predictions for the observable and physical 
properties of galaxies that may be detected by the recently launched Herschel Space 
Observatory in deep fields such as GOOD?>- Herschel. We compare our predictions 
for differential galaxy number counts in the Herschel PACS (100 and 160 jim) and 
SPIRE (250, 350, and 500 jim) bands with available observations. We find very good 
agreement with the counts in the PACS bands, for the overall counts and for galaxies 
binned by redshift at 2: < 2. At 2: > 2 our model underpredicts the number of bright 
galaxies in the PACS bands by a factor of ten. The agreement is much worse for all 
three SPIRE bands, and becomes progressively worse with increasing wavelength. We 
discuss a number of possible reasons for these discrepancies, and hypothesize that the 
effect of blending on the observational fiux estimates is likely to be the dominant issue. 
We note that the PACS number counts are relatively robust to changes in the dust 
emission templates, at least for the three sets of templates that we have tested, while 
the predicted SPIRE number counts are more template dependent especially at low 
redshift. We present quantitative predictions for the relationship between the observed 
PACS 160 and SPIRE 250 /im fiuxes and physical quantities such as halo mass, stellar 
mass, cold gas mass, star formation rate, and total infrared (IR) luminosity, at different 
redshifts. We also present predictions for the radial sizes of i7er5c/ie/-selected disks at 
high redshift {z > 2) and find reasonable agreement with the available observations. 
Finally, we present quantitative predictions for the correlation between PACS 160 
/im fiux and the probability that a galaxy has experienced a recent major or minor 
merger. Although our models predict a strong correlation between these quantities, 
such that more IR-luminous galaxies are more likely to be merger-driven, we find that a 
significant fraction (more than half) of all high redshift IR-luminous galaxies detected 
by Herschel are able to attain their high star formation rates without enhancement 
by a merger. 

Key words: methods: numerical - galaxies: formation - galaxies: evolution - galaxies: 
mergers - galaxies: star formation - galaxies: infrared selected 



1 INTRODUCTION 

One of the most important discoveries from extragalactic ob- 
servations in the mid- and far-infrared has been the identifi- 
cation of luminous and ultra-luminous infrared bright galax- 
ies (LIRGs; LiR > lO^^L© and ULIRGs; Lir > lO^^L©, re- 
spectively). These objects have been studied extensively in 
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the literature (e.g. 


Rieke & Low 1972 Harwit et al. 


1987 


Sanders et al.|1988| 


Sanders, Scoville & Soifer'1991||Condon 


et al.|1991| Auriere et al.|1996| |Duc, Mirabel & Maza 


1997 


Genzel et al. 1998; Lutz et al. 1998; Rigopoulou et al. 
Rowan- Robinson 2000^; ^Genzel et al. 2001^jColina et al. 


1999 
2001 


Colbert et al.||2006| |Dasyra et al.||2006b| |Hernan-Caballero 


et al.||2009| |Magdis et al. 2011), and found to emit 


more 



energy at infrared wavelengths (^ 5 — 500 /xm) than at all 



other wavelengths combined (e.g. [Sanders Mirabel|1996 ). 
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Even though the luminous infrared galaxies are rare objects 
in the local Universe (e.g. Lagache, Puget Dole 2005), 
reasonable assumptions about the lifetime of the infrared 
phase (see, e.g. Farrah et al.|[2003 and references therein) 
suggest that a substantial fraction of all massive galaxies 
pass through a stage of intense infrared emission during their 
lifetime (Sanders Mirabel 1996 , and references therein). 
Consequently, the majority of the most luminous galaxies 
in the Universe emit the bulk of their energy in the far- 
infrared. Furthermore, LIRGs and ULIRGs dominate the 
cosmic star formation rate density at z ^ 1 — 2, account- 
ing for 70 per cent of the star formation activity at these 



epochs (e.g. [Floc'h et al. 2005). This makes the IR an ex- 
tremely interesting wavelength regime to study, especially in 
the context of the cosmic star formation history, and galaxy 
formation and evolution. 

Understanding IR bright galaxies is especially impor- 
tant because light from bright, young blue stars is often 
attenuated by dust (for a review see Calzetti 2001). The 
observed rest-frame ultra-violet (UV) light of a galaxy may 
therefore provide a biased view of the star formation rate 
in the galaxy. For example, the global star formation rate 
required to explain the far-infrared and submillimetre back- 
ground appears to be higher than that inferred from the 
data in the UV-optical (e.g. Hughes et aL][l998 , and refer- 
ences therein). The dust attenuated light from blue stars is 
however not lost, but is re-radiated at IR wavelengths. As 
a result, a large fraction of the radiation from star forma- 
tion is radiated not at UV or optical but at IR rest-frame 
wavelengths. 

The European Space Agency's Herschel Space Obser- 
vatory dPilbratt et"aL]|2010| ) was launched in May of 2009. 
It observes the Universe at infrared wavelengths, from 60 to 
670 /im, and opens up a huge region of new parameter space 
for surveys in area, depth and wavelength. Consequently, a 
primary goal of Herschel is to explore the evolution of ob- 
scured galaxies. With Herschel^ we can finally probe the IR 
light of high-redshift galaxies from 70 to 160 and 250 to 500 
/xm with the Photodetector Array Camera & Spectrometer 
(PACS; |Poglitsch et al.|2010|) and Spectral and Photometric 
Imaging Receiver (SPIRE; Griffin et al.|20ld ), respectively. 
Unfortunately, at such long wavelengths, contamination and 
crowding, rather than the depth of the observations, often 



becomes a limiting factor (e.g. Blain, Ivison k, Smail|| 199*8 



Brisbin et al.||2010| |Lacey et ai.||2010| ). As a result, the in- 
terpretation of IR observations can be less than straightfor- 
ward. This is especially true when probing galaxies at high 
redshift. Despite the potential complications, several deep 
Herschel observations have been performed (e.g. 



Oliver et al. 



2010 ) since the launch of the observatory and more will likely 
follow. 

Early Herschel observations and studies have already 
generated some very interesting results, such as the conclu- 
sions that the majority of the detected high-redshift galaxies 
are large and massive spiral galaxies (e.g. Cava et al.||2010 ) 
with extremely high star formation rates (e.g. Brisbin et al. 
2010). Furthermore, the minimum dark matter halo mass 



of IR bright galaxies was recently constrained to be about 
3 X lO^^M0( [Ainblard et aT]|201l| ), while the average dust 
temperature of H- ATLAS sources was found to be 28 ± 8 K 
( Amblard et al.||2010 ). However, because of the very large 
beam of the Herschel telescope and the difficulty of making 



unique associations between optically detected galaxies and 
the PACS or SPIRE sources, many open questions about 
the physical nature of the populations detected by Herschel 
remain: for example, what are their stellar masses, sizes, 
and morphologies (spheroid or disk dominated)? Are they 
dynamically relaxed galaxies, or interacting or merging sys- 
tems? 

One can turn to theoretical modelling to aid in the in- 
terpretation of these observations, and to make predictions 
about the nature of the objects that Herschel should detect 
at various redshifts. In the IR, all the theoretical models in 
the literature can be divided into two broad classes: phe- 
nomenological models and cosmological models. The phe- 
nomenological models (sometimes called "backwards evolu- 
tion" models) are based on observed luminosity functions, 
and assume simple functional forms to describe the evolution 
of galaxy luminosity functions or SFR with redshift. Tem- 
plate spectral energy densities (SEDs) based on observed 
galaxies are used to transform between different wavelengths 



(e.g. Dale et al.||2001 Dale k Helou 2002i Lagache, Dole k 



Puget 2003 ). In contrast, in semi-analytic cosmological mod- 
els, the evolution of galaxy properties is predicted based 
on the framework of the hierarchical structure formation, 
or Cold Dark Matter (CDM) theory. In these simulations, 
the evolution of structure in the dark matter component 
is characterized via "merger trees" , and simplified prescrip- 
tions are used to model the main physical processes such as 
cooling and accretion of gas, star formation, chemical evolu- 
tion, stellar, supernovae and active galactic nuclei feedback 



(e.g. Croton et al.,2006 i Somerville et al. 2008 ) . Based on the 



resulting predicted star formation and chemical enrichment 
histories, the SED for unattenuated light from stars (which 
dominates the SED shortwards of about 3 /xm) may be calcu- 



& Gunn||1976| Worthey| 


1994| Leitherer k Heckman„1995| 


Bruzual k Charlot||2003 
. B^^^^i : ' 


Maraston||2005| |Conroy, Gunn k\ 



Modelling the impact of dust attenuation, and comput- 
ing the SED of the dust emission, are less straightforward. A 
standard assumption is that all energy absorbed by dust is 
re-radiated. One can then break the problem into two parts: 
computing the amount of light absorbed (and scattered) by 
dust as a function of wavelength, and computing the amount 
of re-radiated light emitted by dust as a function of wave- 
length. One approach is to couple the semi-analytic models 
with a radiative transfer (RT) co de ([Granato et al. ||2000 



Baug h et al.|2005a| |Fontanot et al.||2007| |Lacey et al.||2008 



,2010) to compute the dust attenuation and scattering, and 
a detailed dust model, that assumes a specific dust com- 
position and set of grain properties, to compute the dust 
emission SEDs. Another approach is to use analytic models 
relating the attenuation to the column density of dust (as- 
sumed to be traced by metals in the cold gas) in galactic 
disks (e.g. i Guiderdoni k Rocca-Volmerange|1987 ). To com- 
pute the dust emission SED, one can then use empirical dust 
SED templates based on analytic dust models calibrated to 
observations, or observed galax ies ([Guiderdoni et al.||1998| 
Devriendt, Guiderdoni k Sadat||1999| [Devriendt k Guider-| 
doni||2000t . 

The former approach (full RT+dust model) has the 
clear advantage of providing more detailed and accurate 
predictions of the galaxy SEDs, given the input assump- 
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tions, and certainly provides a better estimate of galaxy-to- 
galaxy scatter in the SED properties. It has the disadvantage 
of being computationally quite expensive, and perhaps not 
merited given that semi-analytic models do not provide de- 
tailed information about the relative geometry of stars and 
dust in galaxies, nor about the composition or properties of 
the dust, which may vary with cosmic time or environment. 



(1999) and described in detail in Somerville et al. (2008), 



Fontanot et al. (2009b) and Fontanot k Somerville (2010) 



have shown that the analytic recipes for predicting dust 
attenuation and the empirical template approach for dust 
emission provide reasonably good agreement, in most cases, 
with statistical quantities such as galaxy luminosity func- 
tions and counts as computed with the full RT approach us- 
ing the GRASIL code ( Silva et al.|1998 ). Moreover, this ana- 
lytic approach to dust attenuation and emission was adopted 
by [Somerville et al.| (2011 hereafter Sll) and |Gilmore et ah] 
(2011), who showed that it provided good agreement with 
observed galaxy luminosity functions for rest-UV to NIR 
wavelengths and redshifts < z < 4. This approach also pro- 
duced fairly good agreement with observational estimates of 
the total IR luminosity function at < ^ < 2, and with the 
observed Extragalactic Background Light. However, these 
models produced poorer agreement at longer wavelengths 
(A ^ 8/im), particularly for luminous higher redshift [z ^ 2) 
galaxies. Thus, the limitations of our approach should be 
kept in mind. Part of the goal of this work is to determine 
where improvements to this kind of analytic approach must 
be made by confronting the models with recent observations 
from Herschel. 

In this work we use the SAM introduced in ISomervillel 



& Primack (1999), with significant updates as presented by 



Somerville et al. (2008), Hopkins et al. (2009), and Sll, and 



undertake a study of galaxies that should be detected by 
current Herschel surveys both locally and at high redshift. 
We focus on the populations that can be detected in the 
relatively deep surveys being carried out with Herschel^ in 
particular the GOODS- if ersc/ie/ project (PI D. Elbaz). The 
redshift range 2 < 2; < 4 is of special interest because 1) it 
spans the peak of the cosmic star formation activity in the 
Universe 2) it probes what are probably the earliest epochs 
for which individual galaxies can be detected with Herschel 
and 3) many open questions remain about the nature of 
galaxy populations, especially luminous IR galaxies, at these 
epochs. 

This paper is organised as follows. In Section [2] we pro- 
vide a brief overview of the semi-analytic galaxy formation 
model used in our study, the model for dust attenuation and 
emission, how we construct our mock light cone, and our 
sample selection criteria. In Section [3] we present our pre- 
dictions for observable properties of Herschel galaxies such 
as counts and luminosity functions. In Section ^ we show 
predictions for the physical properties of Herschel-seXecied 
galaxies at different redshifts. Finally we summarise our re- 
sults and conclude in Section [5] 



2 THE SEMI- ANALYTIC MODEL 

2.1 The Semi- analytic Galaxy Formation Model 

We compute the formation and evolution of galaxies within 
the ACDM cosmology using the semi-analytic galaxy for- 



hereafter SOS. Our model also includes the modified recipe 
for bulge formation and starburst efficiency described in 



Hopkins et al. (2009|. In summary, the SOS SAM includes the 
following physically motivated recipes: 1) growth of struc- 
ture in the dark matter component in a hierarchical cluster- 
ing framework, as characterized by "merger trees"; 2) the 
shock heating and radiative cooling of gas; 3) conversion of 
cold gas into stars according to an empirical "Kennicutt- 
Schmidt" relation; 4) the evolution of stellar populations; 

5) feedback and metal enrichment of the Interstellar (ISM) 
and Intracluster Medium (ICM) from supernovae explosions; 

6) two modes ("quasar" and "radio" mode) of black hole 
growth and feedback from active galactic nuclei; 7) star- 
bursts and morphological transformation by galaxy mergers. 
The SAM can be used to predict physical quantities such as 
stellar and cold gas masses and metallicities of galaxies, their 
star formation and merger histories, and their luminosities 
and sizes. 

The merging histories (or merger trees) of dark 
matter haloes are constructed based on the Extended 
Press- Schechter formalism using the method described in 
Somerville & Kolatt ( 1999|), with improvements described 



in SOS. These merger trees record the growth of dark mat- 
ter haloes via merging and accretion, with each "branch" 
representing a merger of two or more haloes. In this work 
we follow each branch back in time to a minimum progenitor 
mass of lO^M©. This mass scale is sufficient to resolve the 
formation histories of the relatively bright, massive galaxies 
that we study in this work, which occupy halos with masses 
greater than a few xlO"'^"'^M0. 

In the SAM of SOS the star formation occurs in two 
modes, a "quiescent" mode in isolated disks, and a merger- 
driven "starburst" mode. Star formation in isolated disks 
is modelled using the empirical Kennicutt-Schmidt relation 
( |Kennicutt|19S9t , assuming that only gas above a fixed crit- 
ical surface density is eligible to form stars. The efficiency 
and timescale of the merger driven burst mode is a func- 
tion of merger ratio and the gas fractions of the progenitors, 
and is based on the results of hydrodynamic simulations 
( [Robertson et al.||2006| [Hopkins et~aL||2C)09 ). 

In this work we have assumed a standard ACDM uni- 
verse and a Chabrier stellar initial mass function (IMF; 



[Chabrier 2003| ) . We adopt the following cosmological param- 
eter values: = 0.2S, = 0.72, Hq = 70.0, as = O.Sl and 
Us — 0.96, which are consistent with the five year Wilkinson 
Microwave Anisotropy Probe results (Komatsu et al."2 009| 
and was also adopted by Sll. The adopted baryon fraction 
is 0.165S. 

In the following we summarise the main aspects of the 
modelling of the dust attenuation and emission. A more de- 
tailed account can be found in Sll. 



2.1.1 Model for Dust Attenuation 

Our model for dust extinction is based on the approach pro- 



mation model (SAM) introduced in Somerville & Primack 



posed by Guiderdoni & Rocca-Volmerange ( 19S7), combined 
with the model proposed by Chariot Fan[ ( 2000| ) . As in the 
Chariot & Fall model, we consider extinction by two com- 
ponents, one due to the diffuse dust in the disc and another 
associated with the dense 'birth clouds' surrounding young 
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star forming regions. The IZ-band, face-on extinction optical 
depth of the diffuse dust is given by 



rv,o 



^dust,0 ^cold ^cold 



(1) 



where Tdust,o is a free parameter, Zcoid is the metalhcity of 
the cold gas, mcoid is the mass of the cold gas in the disc, 
and rgas is the radius of the cold gas disc, which is assumed 
to be a fixed multiple of the stellar scale length. 

To compute the actual extinction we assign a random 
inclination to each galaxy and use a standard 'slab' model; 
i.e. the extinction in the V-band for a galaxy with inclination 
i is given by: 



-2.51ogi 



exp 



( zixa\ 

\ cos(t) J 



s(i) 



(2) 



Additionally, stars younger than tec are enshrouded in a 
cloud of dust with optical depth tbcv = Mbctv,o, where 
we treat tsc and /xbc as free parameters. Finally, to extend 
the extinction estimate to other wavebands, we assume a 
starburst attenuation curve ( |Calzetti||l997[ (2001 1 for the 
diffuse dust component and a power- law extinction curve 
Ax PC ( A/5500 A ) ^, with n = 0.7, for the birth clouds ( [Char- 
iot k Fall||2000t . Our results are insensitive to minor mod- 
ifications of the underlying extinction curve for the cirrus 
extinction, e.g. adoption of a Galactic or SMC-type extinc- 
tion curve instead of Calzetti. 



2.1.2 Dust Parameters 

There are three free parameters that control the dust at- 
tenuation in our model: the normalization of the face-on 
y-band optical depth rdust,o, the opacity of the birth clouds 
relative to the cirrus component /xbc, and the time that 
newly born stars spend enshrouded in their birth clouds, 
tBC. We first set Tdust,o by matching the normalization of 
the observed relationship between Ldust/Luy vs. bolomet- 
ric luminosity Lboi for nearby galaxies, where Ldust is the 
total luminosity absorbed by dust and re-emitted in the mid- 
to far-IR and Luv is the luminosity in the far-UV 1500 
A). Using a value of Tdust,o = 0.2, Sll found good agreement 
with these observations, and also with the observed optical 
through NIR luminosity functions in the local Universe. We 
also adopt this same value, rdust,o = 0.2. 

The birth cloud parameters /xbc and tBC mainly control 
the attenuation of UV light relative to longer wavelengths. 
Sll show that in the local Universe = 0), the g through 
K-hand luminosity functions are insensitive to the birth 
cloud parameters, while the FUV through i^-bands are quite 
dependent on them. Sll adjusted the parameters to match 
the z = FUV and NUV observed luminosity functions, 
finding good agreement with /xbc = 4.9 and tBC = 2x 10^ yr, 
which we also adopt in our study. However, Sll found that it 
was not possible to reproduce the observed rest-UV and op- 
tical luminosity functions at high redshift with fixed values 
of these parameters. Other studies (e.g. Lo Faro et al.|20()9 



Guo White|200"9| have reached similar conclusions. More- 
over, there is direct observational evidence that galaxies are 
less extinguished at high redshift for a given bolometric lu- 



adopt a simple redshift dependence in all three dust param- 
eters, which is adjusted in order to reproduce the observed 
rest-frame UV and optical luminosity functions out to z ^ 5. 
Our adopted scalings are: rdust,o(^) = T"dust,o(l + and 
both /iBC and tBC scale with z~'^ above z = 1. 



2.1.3 Model for the Dust Emission 

Using the formalism presented above we can compute the to- 
tal fraction of the energy emitted by stars that is absorbed 
by dust, over all wavelengths, for each galaxy in our simu- 
lation. We then assume that all of this absorbed energy is 
re-radiated in the IR (we neglect scattering), and thereby 
compute the total IR luminosity of each galaxy Lir. We 
then make use of dust emission templates to determine the 
SED of the dust emission, based on the hypothesis that the 
shape of the dust SED is well-correlated with Lir. The un- 
derlying physical notion is that the distribution of dust tem- 
peratures is set by the intensity of the local radiation field; 
thus more luminous or actively star forming galaxies should 
have a larger proportion of warm dust, as observations (e.g. 



Sanders & Mirabel 1996) seem to imply. 



There are two basic kinds of approaches for construct- 
ing these sorts of templates. The first is to use a dust model 
along with either numerical or analytic solutions to the 
standard radiative transfer equations to create a library of 
templates, calibrated by comparison with local prototypes. 
This approach was pioneered by [Desert, Boulanger Puget| 
(1990), who posited three main sources of dust emission: 
polycyclic aromatic hydrocarbons (PAHs) , very small grains 
and big grains. In this approach, the detailed size distribu- 
tions are modelled using free parameters, which are cali- 
brated by requiring the model to fit a set of observational 
constraints, such as the extinction or attenuation curves, 
observed IR colours and the IR spectra of local galaxies. 

The second approach is to make direct use of observed 
SEDs for a set of prototype galaxies and to attempt to in- 
terpolate between them (e.g. Chary k, ElbazpOOl Dale & 



Helou 2002 Lagache et al.||2004|). In this work we make 
use of the empirical SED templates recently published by 
Rieke et al. ( [2009 ), hereafter R09. They constructed de- 
tailed SEDs from pubhshed ISO, IRAS and NICMOS data 
as well as previously unpublished IRAC, MIPS and IRS 
observations. They modelled the far infrared SEDs assum- 
ing a single blackbody with wavelength-dependent emissiv- 
ity. The R09 library includes fourteen SEDs covering the 
5.6 X lO^L© < Lir < lO^^L© range. To explore the sen- 
sitivity of our results to the details of the template set, we 
also present results for the Chary Eibaz| ( [2001 ) templates, 
hereafter CEOl, and for the templates recently presented 
by [Chary & Popej (|2010|), hereafter CPU. The latter are 



minosity (Reddy et al. 2010). Following Sll, we therefore 



calibrated to reproduce the SEDs of high redshift galaxies, 
which appear to deviate somewhat from the local templates 
(for a detailed discussion, see CPU). 



2.2 The Mock Light Cone 

We use the SAM described in Section [2T] to generate a mock 
light cone of simulated galaxies. Our simulated light cone is 
100 times the size of the Great Observatories Origins Deep 
Survey (GOODS) on the northern sky, which covers ^ 160 
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arcminutes squared in the sky. We chose to simulate a hght 
cone that covers 100 times larger area than the GOODS-N to 
improve the statistics, as luminous IR galaxies are relatively 
rare. We further limit our light cone of simulated galaxies 
to range from zero to six in redshift space. This limitation 
is arbitrary, however, we will show below that the apparent 
IR flux of simulated galaxies drops steeply as a function of 
redshift, so that it is unlikely that any galaxies beyond this 
redshift could be detected by Herschel. 

Our simulated light cone contains in total 23, 980, 599 
galaxies. The dark matter halo masses range from ^ 2.7 x 10^ 
to ^ 4.5 X lO^^M©. Consequently, our lightcone contains a 
variety of dark matter haloes from light subhaloes to large 
and massive cluster sized haloes. As a result, these simu- 
lated dark matter haloes host a variety of different types of 
galaxies from small dwarf galaxies to giant ellipticals. The 
stellar masses of our simulated galaxies range roughly five 
orders in mae: nitude from lO'^ to - 4.6 x lO^^M©. 

We use the full mock light cone from redshift zero to 
six when comparing the predicted number counts to obser- 



vational constraints (Section 3.1) and for predictions of the 
galaxy luminosity functions (Section 3.2). We also use all 
galaxies in our full mock light cone when making predic- 
tions for physical properties as a function of redshift. How- 
ever, in some cases we select a subsample that contains only 
high-redshift IR bright galaxies. 

2.2.1 Sample of High- Redshift luminous IR Galaxies 

For part of our analysis, we focus on a sample with crite- 
ria that are motivated by the desire to compare with high 
redshift galaxies in the GOODS- ifersc/ie/ survey. For this 
sample, we select all simulated galaxies from our light cone 
that are in the redshift range 2 < z < 4 and above the 
detection limit of either PACS at the 160 /im or SPIRE at 
the 250 /im bands in the GOODS North. For simplicity, we 
assume the limits for both instruments to be ~ 5 mJy (see, 
e.g. Elbaz et al.||2010J , although we note that the actual 
detection limit in the PACS band may be somewhat lower 
(^ 3.5 — 4.5 mJy). These two bands were chosen to represent 
the two different instruments and to be close to the peak of 
the thermal IR emission in our chosen redshift range. No 
other selection criteria, such as stellar mass or optical flux 
cut, were applied, except the ones enforced by the mass res- 
olution. 

Selected in this way, our high-redshift sample comprises 
in total 1154 (3948) galaxies when the selection is done using 
the PACS 160 (SPIRE 250) /im band. Note however that our 
simulated light cone is 100 times larger than the GOODS-N 
— hence one may infer that we would expect ^ 10 PACS 
160 /im and ^ 40 SPIRE 250 /xm detected galaxies in a 
field with the size and depth of GOODS-N. However, cosmic 
variance is probably large for these objects, and therefore 
this number is uncertain. We return to this matter in the 
following Sections. 

The PACS 160 /xm selected high-redshift sample con- 
sists of luminous IR galaxies with bolometric IR luminosity 
ranging from ^ 1 x 10^^ to ^ 2 x lO^^L©. The average bolo- 
metric luminosity Ljj^ - 2.46 x lO^^L©. The SPIRE 250 /im 
selected high-redshift sample consists of IR bright galaxies 
with the bolometric IR luminosity ranging from ^ 5 x 10^^ 
to ^ 2 X lO^'^L©, while the average bolometric IR luminos- 



ity Ljp^ ^ 1.75 X lO^^L©. We can therefore readily conclude 
that all of the galaxies in the PACS selected sample are 
ULIRGs, while in the SPIRE selected sample they are at 
least LIRGs, while most, i.e. more than 90 per cent, can 
be defined as ULIRGs. The simple single PACS 160 /xm or 
SPIRE 250 /im band selection is a good indicator for the 
total IR luminosity because in our selected redshift range 
we are always near the peak of the IR flux. Consequently, 
our sample of high-redshift galaxies is well suited to study 
the physical properties of high redshift (U) LIRGs that are 
potentially observable in fields such as the GOODS North 
and South. 

Figure[l]shows the PACS 160 /im flux, 5*160 , in mJy as a 
function of redshift z for all simulated galaxies as well as the 
co-moving number density of simulated galaxies as a func- 
tion of PACS 160 /im flux and redshift. The horizontal green 
dashed line in the scatter plot is drawn at 5 mJy, which is 
roughly the same as the detection limit in the GOODS-N. 
Note, however, that the confusion limit in the GOODS-N 
and in comparable fields can be higher than the detection 
limit, especially at longer wavelengths. For example, a con- 
fusion hmit - 19.1 ± 0.6 mJy at the 250 /xm SPIRE band 



has been quoted in literature (Nguyen et al. 2010). Obvi- 



ously, observations that are of the same depth, but on less 
crowded fields could reach down to the mJy level. Moreover, 
observations in gravitational lensed fields can reach well be- 
low the derived confusion limit, as has already been shown 
by |Altieri et al. (2010), thus providing a motivation for our 
sample selection and the usage of 5 mJy. We note that the 
Fig. [l] would be quantitively similar if we used the SPIRE 
250 /im band instead of the PACS 160 /im flux, although the 
SPIRE-selected sample extends to slightly higher redshift (a 
few galaxies are predicted above 5 mJy up to 2; ^ 4.0). 

Figure [l] shows that the number density of luminous 
IR galaxies drops quickly above redshift three, hence, ob- 
servations comparable to those in GOODS-N are unlikely to 
detect galaxies above redshift four. However, as noted above, 
these galaxies are likely to be strongly clustered, with a large 
field-to- field variance, implying that the maximum redshift 
may also vary from field to field. Unfortunately, due to the 
fact that we adopt the Extended Press- Schechter formalism 
to generate our merger trees, we do not have information 
on the spatial locations of our galaxies and therefore we 
cannot properly quantify the effects of cosmic variance. We 
defer this to a future study, where we plan to use lightcones 
extracted from A/'-body simulations and will therefore have 
information on the clustering properties of our galaxies. 



3 OBSERVABLE PROPERTIES 

3.1 Number Counts and Redshift Distributions 

The most basic statistic describing a galaxy population is 
the number counts, i.e. the number density on the sky of 
galaxies as a function of observed flux. The number counts 
at far-infrared and sub-mm wavelengths are well known to 



exhibit strong evolution ( Oliver et al.||2010 and references 
therein). Unfortunately, the Herschel beam is broad com- 
pared to the number density of sources, and thus the maps 
are often confused. This confusion means that care has to be 
taken when estimating number counts from Herschel data: 
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Figure 1. PACS 160 fim flux, Siqq, in mJy as a function of redshift z for aU simulated galaxies. The histograms show projections of the 
scatter plot. The solid bars in the upper histogram show co-moving number density (^^) of all galaxies, while the hatched bars show that 
of simulated galaxies with Siqq > 5 mJy. The right-hand side histogram shows co-moving number density (^) as a function of PACS 160 
fj,m flux. The horizontal green dashed line is at 5 mJy, which corresponds roughly to the detection limit in PACS 160 fim observations in 
the GOODS-N. 



observations must in general be corrected for flux boosting 



and incompleteness (see, e.g., [Clements et al. 2010 Oliver 
et al.|[20To and references therein). For example, the faint 



flux densities may be overestimated due to the classical flux 
boosting effect (e.g. |Bethermin et al.|2010 ). In the following 
sections we therefore compare our model to existing obser- 
vations and make predictions for higher redshifts and lower 
flux levels than what the current observations have probed. 



3.1.1 Number Counts in the PACS and SPIRE Bands 

Figures |2] and [3] show diflFerential galaxy number counts for 
the PACS 100 and 160 fim. bands derived from our simu- 
lated light cone after being normalised to a Euclidean slope 



(dN/dS oc S- 



In order to investigate the sensitivity of 



our results to the dust emission templates, we have shown 
the predictions using three different sets of templates. The 
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10° IQi 10^ 

S'loo [mJy] 

Figure 2. Predicted galaxy differential number counts (circles) 
compared to observational estimates of |Berta et al.| ( |2010| ) (shown 
in cyan) and |Altieri et al^ ( |2010| ) (shown in magenta) in the PACS 
100 /im band at different redshift bins. The |Berta et al.| (|2010[) 
observations are based on direct detections, while the |Altieri et al.| 
( |201 0) number counts take advantage of lensing caused by the 
cluster Abell 2218. The shaded regions correspond to 3cr Poisson 
errors in the model, while the observational errors are Icr limits. 
The different colours refer to different dust emission template 
SEDs that have been used in the semi-analytic model (red: [Chary 




I&: Elbaz, ( ,2001^ ; blue: , Chary Pope, ( ,2010, ); green: [Rieke et al. 
(|2009|> ). 



10° IQi 10^ 

5*160 [mJy] 

Figure 3. Predicted galaxy differential number counts (circles) 
compared to observational estimates of ,Berta et al.| ( |2010^ (shown 
in cyan) and |Altieri et al^ ( |2010| > (shown in magenta) in the PACS 
160 /im band at different redshift bins. The [Berta et al.| (|2010|> 
observations are based on direct detections, while [Altieri et al.| 
(2010) number counts take advantage of lensing caused by Abell 
2218. The shaded region corresponds to 3cr Poisson errors in the 
model, while the observational errors are Icr limits. The different 
colours refer to different dust emission template SEDs that have 
been used in the semi-analytic model (red: Chary & Elb az| ( |^001| >; 
blue: |Chary fc Pope, ( ,2010, ); green: |Rieke et al.| ( |2009| ) ). 



template sets of [Rieke et al. ( 2009 ) and Chary & Elbaz 
( 2001 ) are based on nearby galaxies. The templates oflChary 



& Pope (2010) incorporate data from distant galaxies, and in 



general predict higher fluxes longwards of the thermal dust 
peak for IR luminous galaxies, corresponding to colder dust 
temperatures. Overplotted are the galaxy differential num- 
ber counts from the early Herschel observations of [Berta[ 
et al.[ ([2010[) and [Altieri et al.[ ([2010|). The number counts 



of [Berta et al.[ ([2010|), which were derived from GOODS- 
N data, cover the flux range 3 — 50 mJy at 100 /im and 
5.5 — 72 mJy at 160 /im, while the number counts of |AltierTl 
et al. (2010) cover the flux range ^1 — 35 and ^ 1.5 — 58 



mJy at 100 and 160 /im, respectively. The general agreement 
between the galaxy differential number counts of simulated 
galaxies and observations is excellent at low {z < 2) redshift. 
However, in the highest redshift bin (2 < z < 5) the agree- 
ment is worse. Even though the observational constraints are 
weaker — only two (four) of the observational constraints 
at the PACS 100 (160) /im band are not upper limits — this 
disagreement requires further discussion. We return to this 
issue later. 

A detailed inspection of the differential number counts 
of the simulated galaxies shows a monotonic rise in the two 
lowest redshift bins till ^ 20 mJy. This is consistent with a 
non-evolving population of galaxies. Based on their observa- 
tions ,Berta et al. (20 10 J argue for a drop at ^ 20 m Jy, even 
though the errors in observations are relatively large. The 
number counts of the simulated galaxies provide some sup- 
port for such a claim at the 100 /im band, however, at the 
160 /im band our model does not show any significant drop. 
For the PACS 100 /im band our model predicts a reasonably 
monotonic rise all the way to ^ 500 mJy. Even so, a small 



(but statistically significant) drop is visible at ^ 60 mJy. 
Quantitatively the PACS 160 /im band's behaviour is very 
similar, however, the increase in the Euclidean normalised 
number counts is shallower after ^ 20 mJy and seems to 
reach a plateau. The Euclidean normalised counts seem to 
rise again at ^ 200 mJy, especially at the PACS 100 /im 
band. 

The redshift evolution of the galaxy differential number 
counts of the simulated galaxies in both PACS 100 and 160 
/im bands is rather modest at z ^ 1. Even up to redshift 
two the number of ~ 10 mJy IR galaxies seems to evolve 
only slowly, in contrast to the number of fainter galaxies, 
which seems to grow as a function of redshift, causing the 
distribution of the number counts to flatten and eventually 
to turn over when probing even earlier cosmic times. At 
z ^ 2 our model predicts that the normalised number counts 
of the faintest {Sieo ^0.1 mJy) galaxies are at the same level 
as that of the brightest (S'leo 100 mJy). 

Figures ^ [5] and [6] show Euclidean normalised differ- 
ential galaxy number counts derived from our simulated 
light cone in the 250, 350, and 500 /im SPIRE bands, re- 
spectively. We show the observational results from the early 
Herschel observations of Glenn et al. ( 2010, ) and , Clements] 



et al. (2010) for comparison. At 250 /im, our model predicts 
significantly fewer galaxies at the 5 — 30 mJy level than are 
observed. This discrepancy is even larger in the SPIRE 350 
and 500 /im bands. The agreement seems better at slightly 
higher fluxes (~ 100 mJy), but then the models overpredict 
sources at the highest flux levels (^ 1 Jy) in all three bands. 
Note however that the number of simulated galaxies with 
S ^ 400 mJy is very small, as indicated by the large error- 
bars. The actual field-to-field variation will be even larger, as 



© 2011 RAS, MNRAS 000,[lpOl 



8 Sami-Matias Niemi et al 



Total 




0.0 < z < 1.0 



10° IQi 102 103 10° IQi 102 10^ 

5*250 [mJy] ^250 [mJy] 

Figure 4. Predicted galaxy differential number counts (circles) 
compared to observational estimates of |Glenn et al-l i pOlOl ) (shown 
in cyan) and Clements et al. ( 2 Q10| (shown in magenta) in the 
SPIRE 250 /im band in different redshift bins. The number counts 
oflClements et al.|(|2010| are direct detections, while Glenn et al. 
( |2010| derived their counts by modelling the P(D) distribution, 
rather than identifying individual galaxies. The shaded regions 
correspond to 3a Poisson errors in the models, while the obser- 
vational errors are la limits. Green, red, and blue colours show 
the predictions using different dust emission templates, as in the 
previous figures. 




10° 10^ 102 103 10° 10^ 102 103 

^350 [mJy] 5350 [mJy] 

Figure 5. Predicted galaxy differential number counts (circles) 
compared to observational estimates of |Glenii et al.| ( |201Q] (shown 
in cyan) and [Clements et al.] ( |2010| (shown in magenta) in the 
SPIRE 350 /im band at different redshift bins. The number counts 
of Clements et all(|2010|) are direct detections, while [Glenn et ah] 
(2010) derived their counts by modelling the P(D) distribution, 
rather than identifying individual galaxies. The shaded regions 
corresponds to 3a Poisson errors in the models, while the obser- 
vational errors are la limits. Green, red, and blue colours show 
the predictions using different dust emission templates, as in the 
previous figures. 



our model errorbars do not include the effects of large-scale 
clustering. 

The galaxy differential number count estimates of [Glenn] 



et al. (2010) are based on the pixel brightness distribution 



P{D) method, in which no direct detections are needed. In 
this method the distribution of pixel brightnesses are used 
to make statistical interferences about the slope and ampli- 
tude of the source counts. Therefore, the P{D) method can 
be used to derive number counts below the confusion limit 



(see Glenn et al. (2010) for details). However, we note that 



Glenn et al. (2010 ) emphasise that their values are model fits 



and thus effectively integral constraints over some region 
surrounding each flux density. However, as the differences 
between their multiply-broken power-law and spline mod- 
els are relatively small, we choose to show only their results 
based on the multiply-broken power-law model (Table 2 of 
Glenn et al.[[2010|. This provides a reasonable compromise 
given the fact that we have added their la statistical errors 
and their estimated systematic uncertainties in quadrature, 
after which both of their model fits are roughly within the 
errors shown. 

Glenn et al.| (|2010| note that they clearly detect a break 



in the number counts at around 20 mJy in all the SPIRE 
bands at high significance. Our model also shows a break, 
however, it is at a significantly lower flux (^ 2 mJy), except 
for the SPIRE 250 /xm channel when adopting the [Ghary k\ 
Pope (2010) templates. Glenn et al. (2010) also note that 



there is possible weak (~ la) evidence for a 'bump' in the 
differential counts around 400 mJy at the 250 /xm band. 
The location of such a bump in our models is quite depen- 
dent on the adopted dust emission templates, but our model 
shows an upturn at high fluxes independent of the templates 
adopted. However, at the 350 /xm band all templates show a 




10° IQi 102 10^ 10° 10^ 102 10^ 

5*500 [mJy] 55oo [mJy] 

Figure 6. Predicted galaxy differential number counts (black cir- 
cles) compared to observational estimates of [Glenn et al.[ ( |2010[ ) 
(shown in cyan) and [Clements et al.] ( [2010[ > (shown in magenta) 
in the SPIRE 500 /xm band at different redshift bins. The number 
counts of C lements et al.[ ( [2010| ) are direct detections, while Glenn] 
et al. (20101 derived their counts by modelling the P{D) distri- 
bution, rather than identifying individual galaxies. The shaded 
regions correspond to 3a Poisson errors in the models, while the 
observational errors are la limits. Green, red, and blue colours 
show the predictions using different dust emission templates, as 
in the previous figures. 
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bump around 300 mJy, although it is unclear whether this 
is part of the upturn or a real feature given the relatively 
large uncertainties at high fluxes. 

We are not aware of any published results for SPIRE 
counts broken into redshift bins, however we show our theo- 
retical predictions for comparison with future observational 
results. The predicted redshift evolution of the number 
counts in all SPIRE bands is qualitatively similar to that 
of the PACS bands. We note however that in the SPIRE 
bands the redshift evolution of LIRGs is more modest. Even 
in the highest redshift bin, the Euclidean normalised galaxy 
differential number counts of the brightest galaxies are only 
a factor of a few lower than in the lower redshift bins. 

In summary, our models produce very good agreement 
for the overall counts in the PACS 100 and 160 /xm bands, 
but this agreement rapidly worsens as we move to longer 
wavelengths. The overall agreement in the SPIRE bands is 
quite poor. Moreover, although our model predictions are 
also in good agreement with the galaxy number counts split 
by redshift for z < 2, they do not produce enough IR- 
luminous galaxies at z > 2. Although the magnitude of the 
discrepancy between our theoretical predictions and the ob- 
servational results is disturbing, the interpretation is also 
somewhat complex. In what follows we therefore discuss 
potential theoretical and observational complications that 
could lead to the noted discrepancies in turn. 

We first emphasise that these semi-analytic models have 



been shown elsewhere (see Sll; Fontanot et al. j2009a, and 
references therein) to reproduce reasonably well observed 
luminosity functions from the rest UV to near-IR, as well 
as stellar mass functions and star formation rate functions, 
from z ~ — 5, although the agreement does get worse to- 
wards higher redshifts. Moreover, the good agreement in the 
PACS bands up to z 2 suggests that many of the basic 
ingredients of the model must be reasonable. However, a 
possible explanation for some of these discrepancies is that 
perhaps the adopted dust emission templates are not valid 
at earlier cosmic times. Indeed, there is observational evi- 
dence that the correlation between total IR luminosity and 
dust temperature might be different in high redshift galax- 
ies ([Amblard et al.||20TQl |Elbaz et al.||2010| [Hwang et aT] 
2010). We attempted to investigate this by testing three dif- 
ferent sets of dust emission templates (see also the discus- 
sion in CPU), as described at the beginning of this section. 
The results of adopting the different templates are complex, 
with some templates producing better results at some wave- 
lengths and redshifts, and others performing better in other 
regimes. The predictions with the three templates are similar 
in the PACS bands for the overall counts and at low redshifts 
{z < 1). The CPU templates predict slightly fewer galaxies 
in the highest redshift bin (2 < z < 5), which worsens the 
(already poor) agreement with the observations. However, 
the CPU templates produce significantly better agreement 
for the counts in the SPIRE bands at intermediate fluxes 
(between about 3 and 100 mJy). Thus, we conclude that the 
results are sensitive to the choice of dust emission templates, 
but none of the template sets that we investigated can repro- 
duce the observations at all wavelengths. Another potential 
cause for the noted discrepancy may relate to our simplified 
modeling of the dust absorption and emission. We tuned our 
model to reproduce the observed ratio of unattenuated to re- 
radiated light (UV to FIR) for average galaxies at z = 0. 



However, galaxies at high redshift may have more complex 
geometries, leading to greater scatter or even a breakdown 
in this approach. As already noted, scatter or evolution in 
the dust emission SEDs could also be important. As noted 
by CPU, it is not clear whether the observations they used 
span the entire range of underlying dust temperatures and 
associated far-infrared colours. 



On the other hand, one likely culprit in the Herschel 
observations is the impact of flux boosting and blending, 
which affects the counts more severely at longer wavelengths 
and may also have a larger effect at high redshifts, where 
galaxies are strongly clustered and there is a higher proba- 
bility that the Herschel beam intersects lower redshift ob- 
jects along the line of sight. Indeed, we have visually checked 
the Hubble Space Telescope (HST) images at the position 
of 250 /xm detected Herschel sources in GOODS-N, and it 
is clear that in many, if not most cases, there are multi- 
ple galaxies within the PACS or SPIRE beam. For longer 
wavelengths such as the SPIRE 500 /im the complications 
are even worse because of the full-width-at-half maximum 
of the beam doubles (from ^ 18 at 250 /xm to ^ 36 arc sec- 
onds). Moreover, according to the simulations carried out 
by Rigby et al. (2011), more than half (^ 57 per cent) of 
sources detected at > 5cr at 500 /xm show a flux boosting by 
a factor > 1.5 and more than every fourth (^ 27 per cent) 
by a factor > 2. This is clearly a very serious issue when 
comparing theoretical predictions with these observations. 



Other observational complications include cosmic vari- 
ance. Particularly the bright counts may be compromised 
by field-to- field variance, as the areas covered are still not 
large. Moreover, at long wavelengths, recent results imply 
that many of the bright sources might be lens ed (e.g. |Ne- 
grello et al.1[20Tol [Lima, Jain Devhn||20Tol [Lima et al." 
2010| ). Finally, the contribution to the dust heating from 
obscured active galactic nuclei (AGN) is uncertain (e.g. |AL 



maini, Lawren ce Boyle P1999!* Granato, Danese Frances- 
chini 1997; Farr ah et a l. 2007; Saj ina et al.||2008| ), although 
this is unlikely to contribute significantly particularly at the 
longer wavelengths. Because we cannot rule out these effects, 
for the remainder of this work, we therefore caution that the 
predicted IR fluxes that we quote should not be interpreted 
as corresponding in a one-to-one fashion with the measured 
fluxes for individual sources that can currently be detected 
in a Herschel image, but rather with the idealised flux that 
one would obtain if flux boosting, blending, and cosmic vari- 
ance could be overcome. This situation is less than ideal, 
but we feel that it can provide a first step towards provid- 
ing useful predictions. In future work, we plan to pursue 
a more rigorous comparison, in which we use the theoreti- 
cal simulations to produce mock Herschel images and treat 
them in the same way as the real observations. Moreover, 
the GOODS- Herschel team and others are also working on 
a better understanding of the impact of blending on their 
observations, which will eventually clarify the interpretation 
of our predictions. In addition, our predictions in this form 
can be compared directly with other theoretical predictions 
in the literature, which similarly do not address the effects 
of blending. 
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3.2 Galaxy Luminosity Functions 

While simple number counts are the most basic statistic of 
galaxy populations, their predictive power is limited. A step 
further towards the intrinsic properties of a galaxy popula- 
tion can be taken by studying the galaxy luminosity func- 
tions (LFs). However, even more knowledge of galaxy evo- 
lution can be gathered when the evolution of the rest-frame 
LFs are studied as a function of redshift. 

Figure |7| shows the rest-frame galaxy LFs in the PACS 
100 and 160 /xm and SPIRE 250 and 350 /im bands. The rest- 
frame galaxy LFs are shown in five different redshift bins: 
< z < 0.3, 1.2<z< 1.6, 2.0 < z < 2.4, 2A < z < 4.0, and 
4.8 < z < 5.1. This enables the study of the redshift evo- 
lution of the galaxy IR LFs between z ^ 5 and z ^ 0. Un- 
fortunately the observational constraints in the rest-frame 
LFs of Herschel bands are still limited, especially at higher 
redshifts, thus, the comparison between the theoretical pre- 
dictions and observations is less than robust. However, to 
provide an initial comparison we show the LFs of Lapi et al.| 
( |2011| , which have been derived from early Herschel data. 
The redshift bins of the theoretical LFs have been matched 



to the ones presented in Lapi et al. ( 2011 ) to allow easy com- 



parison. However, the sample of Lapi et al. (2011) is limited 



to only bright galaxies with S250 > 35 mJy Moreover, even 
though they estimate that the fraction of false identifica- 
tions is as low as 6 per cent (in areas with good SDSS 
completeness) , they warn that some truly high-z sources can 
be missed by their study leading to incompleteness as high 
as 20 per cent. Such complications will obviously make the 
comparison less than robust. 

The upper panels of Figure |7| namely a and b, show the 
predicted galaxy LFs in the PACS 100 and 160 /xm bands. 
Qualitatively these two panels are rather similar, i.e., there 
are only small differences in the LFs of the rest-frame 100 
and 160 /xm bands. We do note however that the 160 /xm 
band does not extend to as high luminosities as the 100 /xm 
band. This is true at all redshifts. The rest-frame LF of the 
lowest redshift bin is the lowest curve, showing that the red- 
shift evolution in the PACS bands is strong. The two panels 
show that the number density of the PACS bright galaxies is 
highest at redshifts z ^ 2, while the highest number density 
of the faint galaxies is found at z ^ 3. Both PACS bands 
show that the faint end of the LFs evolves hardly at all be- 
tween z ^ 5 and z ^ 2. Even though we do not see much 
evolution in the faint end of the PACS LFs, a clear differ- 
ence exists between the local and higher redshift galaxies. 
Our model predicts that the number density of luminous in- 
frared galaxies peaks around z ~ 2 — 3, in agreement with 
other predictions (e.g. Swinbank et al. 2008| Lacey et alj 



2010). Furthermore, the total number of IR bright galaxies 
is highest at 2; ~ 2. Such results are easy to understand if the 
IR rest-frame flux correlates with dust obscured star forma- 
tion rate, as the cosmic star formation rate density is known 



to peak around the same cosmic time (e.g. [Bouwens et al. 
[2011 ) . Our model also predicts that the number density of 
the most luminous IR galaxies (i.e. ULIRGs) increases by a 
factor of ^ 5 between redshifts five and three, while their 
number density is found to peak at z ^ 2. The evolution of 
the PACS LFs also shows that the number density of lumi- 
nous IR galaxies drops by a factor of ^ 10 between redshifts 
two and zero. These results emphasise the strong evolution 



in the number density of (U)LIRGs as a function of cosmic 
time. 

In panel a of F igure [7| we a lso compare our theoretical 



LFs to the LFs of Lapi et al. (2011). Qualitatively simi- 



lar redshift evolution is seen in both cases. Our prediction 
for the lowest redshift observed LF (1.2 < z < 1.6) is in 
reasonable argeement with the observations, with a modest 
over-prediction of the highest luminosity sources. Towards 
higher redshift, and especially in the highest redshift bin, the 
disagreement grows rather large in sense that our model pre- 
dicts fewer galaxies with high luminosities than are reported 
in the observations. Similar effects as discussed in the case 
of number counts are likely to contribute. The most strik- 
ing difference, however, between our model and the LFs of 



Lapi et al. (2011 ) is the location of the knee of the LFs. Our 



predicted LFs turn over at about two orders of magnitude 
higher number densities. One potential explanations is pro- 



vided by Lapi et al. (2011) who argue that the flattening of 
their LFs at the lowest luminosities may be, at least in part, 
due to the overestimate of the accessible volume yielded by 
the l/Vmax estimator for objects near to the flux limit. 

The lower panels of Figure |7| namely c and d, show the 
predicted rest-frame galaxy LFs in the SPIRE 250 and 350 
fim bands. The redshift evolution of the rest-frame SPIRE 
band LFs is similar to that of the PACS bands. Again, the 
number densities of LIRGs and ULIRGs, i.e., the galaxies at 
the bright end of the LF, are found to peak at z ^ 2 — 3. In 
agreement with the PACS bands, the highest number densi- 
ties are found at redshifts two and three, while the number 
densities drop by a factor of ^ 11 to redshift zero. Thus, 
the LFs of the rest-frame SPIRE bands show qualitatively 
similar evolution as the LFs of the PACS bands. However, 
we note that the bright end of the LFs in the SPIRE bands 
appears to drop more steeply than in the PACS bands. 

In panel c of Figure [7|we also compare our theoretical 
LFs to LFs of |Lapi et ari ( |2011| ). As in case of the PACS 100 
/xm the faint end of the lowest redshift LFs (1.2 < z < 1.6) is 
in good argeement with our prediction, while in the brightest 
end our model predicts more galaxies than are observed. 
Towards higher redshift, and especially in case of the highest 
redshift bin, the disagreement grows rather large. Again, 
similar effects as discussed earlier are likely to play a role. 

Infrared LFs have also been derived from other obser- 
vations (e.g . |Kim Sanders||1998| [Floc'h et ar]|2005| [Goto 
et al.pOlO l) and early Herschel data. |Dye et al.| ( |2010t de- 
rived SPIRE 250 fim band LFs out to z = 0.5. They find 
that the LF exhibits significant evolution out to z — 0.5 
and that at a given luminosity, the comoving space den- 
sity increases steadily with redshift. The evolution noted in 
Dye et al. (2010) is qualitatively in agreement with our LFs 



(see Sll), although, their evolution is rather strong given 
the rather narrow redshift range. However, strong evolu- 
tion, particularly at the bright end, is in agreement with 
our LFs and also noted by others (e.g. |Eales et al.||2009) 



Also Vaccari et al. (2010) derived LFs for local galaxies in 



all three SPIRE channels. Their LFs are also in good qualita- 
tive agreement with our findings for local galaxies. However, 
our lightcones have very small volumes at low redshift, as 
they were intended for high redshift comparisons, so a more 
quantitative comparison is beyond the scope of this study. 
A comparison with the z — 0.5 SPIRE 250 jim. LF is shown 
in Sll. 



© 2011 RAS, MNRAS QQQ,[TpOl 



Physical properties of Herschel selected galaxies 11 



I 

X 

CD 

CO 

I 

O 



O 



I 

X 

CO 
I 

o 



O 




9.5 10.0 10.5 11.0 11.5 12.0 12.5 13.09.5 

logio(^ [L©]) 



10.0 10.5 11.0 11.5 12.0 12.5 13.0 

logic [L©]) 



Figure 7. Rest-frame galaxy luminosity functions in the PACS 100 and 160 /im bands and SPIRE 250 and 350 /xm bands at five different 
redshifts. Different redshifts are shown in different colours and are noted in the legend. The observed luminosity functions (in a and c) 
marked with symbols are from [Lapi et al.] ( |2011| . 



4 PHYSICAL PROPERTIES OF IR-LUMINOUS 
GALAXIES OVER COSMIC TIME 

Having explored some of the statistical properties of Her- 
schel detected galaxies, in this Section we present predic- 
tions for some of the physical properties of these objects 
at different redshifts. We explore the correlations between 
160 and 250 /xm flux and many different physical properties 
of our model galaxies. Here we show some of the strongest 
and most interesting correlations that we identified, which 
include stellar mass, dark matter halo mass, cold gas mass, 
star formation rate, and total IR or bolometric luminosity. 
The correlations with the other Herschel bands are similar 
and all qualitative conclusions would remain the same, so 
we do not show them. 

Figure [S] shows the predicted median stellar mass as a 
function of 160 and 250 /im flux, along with the 16th and 
84th percentiles, for different redshifts. Our model predicts 
that all galaxies with 5*160 > 5 or 5250 > 5 mJy have stel- 
lar masses in the range ^ 10^ — 10^^ M©. Furthermore, 
we predict that high redshift {z > 2) galaxies that can be 
detected individually in deep Herschel surveys such as the 
GOODS are quite massive, > lO^^-^M©. This is quite 



different from the results of Lacey et al. (2010), who pre- 



dict much smaller stellar masses for galaxies at a given FIR 
flux — at z ^ 2, even their most luminous galaxies are less 
massive than ^ IO^^Mq. This difference probably arises 



mainly from the fact that Lacey et al. (2010) adopted a 



top-heavy IMF during starbursts. Therefore there are fewer 
long-lived low mass stars in galaxies at high redshift. 

Not surprisingly, as stellar mass is know to be fairly 
well correlated with dark matter halo mass, we find sim- 
ilar results for the dark matter halo masses. At the low- 
est redshifts {z < 0.5), luminous IR galaxies (with 5i6o 
or 5250 > 5 mJy) can be found in dark matter haloes as 
light as logio(Mdm/M0) 11.1 (see Fig. [9|. However, the 
bulk of simulated IR bright galaxies with 5i6o or 5250 > 5 
mJy, especially at higher redshifts, were found to reside 
in relatively massive dark matter haloes. To be more pre- 
cise, our model predicts that the masses of the dark matter 
haloes of IR bright galaxies can cover a broad range from 
log^Q(Mdm/M0) ^ 11.5 to 13.5, while the bulk of luminous 
IR galaxies can be found to reside in rather typical dark 
matter haloes with masses Iog2o(^dm/M0) ^ 12.5. 

These predictions are in good agreement with the re- 



© 2011 RAS, MNRAS 000,[lpOl 



12 Sami-Matias Niemi et al. 



cent results of Amblard et al. (2011), which were derived 



from Herschel data. In this paper the authors find an excess 
clustering over the linear prediction at arcminute angular 
scales in the power spectrum of brightness fluctuations at 



250, 350, and 500 /xm. From this excess [Amblard et al. (2011 ) 
infer that IR bright galaxies are located in dark matter halos 
with a minimum mass of log]^Q(Mmin/M0) = 11.51q"2 ^50 
/xm. When the authors average over the three wavelengths 
the minimum halo mass for their galaxies is at the level of 
3 X 10^ ^Mq. All our IR bright galaxies are found to reside in 
dark matter haloes more massive than this lower limit, ex- 
cept a handful of local galaxies, which reside in haloes only 
slightly less massive (1 x lO^^M©) than the minimum mass 



of Amblard et al. (2011) 



Our simulated galaxies show evidence for a weak cor- 
relation between the dark matter halo mass and the PACS 
160 and SPIRE 250 /xm fluxes (Fig.|9|. Statistically, galax- 
ies residing in more massive dark matter haloes have higher 
median IR flux as was also the case with stellar mass (Fig. 
|8|. However, the trend between the dark matter halo mass 
and 160 or 250 /xm flux is weaker than in case of stellar 
mass as indicated by the larger range between the 16th and 
84th percentiles in Fig |9] Moreover, for the high redshift 
galaxies (2 < z < 4) with either 5*160 or 5250 > 5 mJy 
no statistically significant correlation can be inferred. In- 
stead, galaxies with the highest fluxes are again found to 
reside in dark matter haloes with a broad range of masses 
from log^Q(Mdm/M0) ^ 12 to 13.5. It is however interest- 
ing to note that luminous IR galaxies can also be found 
in cluster-mass dark matter haloes. This implies that some 
luminous IR galaxies at high redshift may be the progeni- 
tors of present-day brightest cluster galaxies (BCGs). More- 
over, this also implies that the IR bright galaxies at high 
redshifts may provide a signpost to identify proto-clusters 
during the epoch of their formation. Note, however, that ac- 
cording to our model, the PACS 160 and SPIRE 250 /xm 
fluxes are almost independent of the dark matter halo mass 
for high-redshift galaxies that would be detected in Herschel 
observations comparable to fields such as the GOODS. Con- 
sequently, the IR flux at high redshift must be governed pri- 
marily by other quantities than the mass of the dark matter 
halo. 

Figure [To] shows the predicted median mass in cold gas 
as a function of the PACS 160 and SPIRE 250 /xm flux for 
galaxies selected at different redshifts. Not surprisingly, our 
model predicts that at high redshift the luminous IR galaxies 
directly detectable in Herschel surveys will be the most gas 
rich. Our simulation shows that the galaxies found in these 
surveys should typically be gas rich with cold gas masses 
ranging from - 10^ up to - 4 X lO^^M©. The bulk of the 
high redshift IR bright galaxies were found to contain a cold 
gas mass of - 6 X 1O^°M0. Interestingly, as we noted earlier, 
the bulk of the high redshift IR bright galaxies have stellar 
masses 6 x lO^^M©. Consequently, the high-redshift 

luminous IR galaxies should be extremely gas rich, with the 
mass in cold gas comparable to or even higher than the mass 
in stars for high-redshift luminous IR galaxies. This is a clear 
prediction for future observatories such as ALMA, which 
should be able to detect CO emission from the molecular 
gas in these galaxies. 

Figure [TT] shows the predicted median star formation 
rate as a function of the PACS 160 and SPIRE 250 /xm 



flux at different redshifts. It is not surprising that this is 
the strongest correlation seen so far, as the far-IR flux is 
generally assumed to be a good proxy for the star forma- 
tion rate, and at some level this correlation is artificially 
tight because of our assumed one-to-one correspondence be- 
tween total IR luminosity (light reprocessed by dust) and 
dust emission SED. However, this plot provides a conve- 
nient guide to what SFR one expects to be able to power 
galaxies at a given 160 or 250 /im flux at different redshifts, 
assuming a standard set of templates, as well as an indica- 
tion of how much scatter can arise from heating of dust by 
old stars. 

Figure [TT] shows that at z > 2 the galaxies directly 
detectable in surveys such as the GOODS-N should have 
median star formation rates > 3OM0yr~^ according to our 
model prediction. We find that the average star formation 
rate for luminous IR galaxies at a high redshift (2 < 2; < 4) 
is ^ 155M0yr~^. This is significantly higher than for all 
galaxies (with log;Lo(^dm) > 9) in the same redshift range, 
for which our model predicts a mean star formation rate of 
merely ^ 2M0yr~^. The highest star formation rate pre- 
dicted by our model for a high-redshift galaxy is as high as 
^ 78OOM0yr~^. In the next Section, we explore the physi- 
cal mechanisms that are driving high SFRs in high redshift 
galaxies. Clearly, such high SFR cannot persist for very long 
periods of time! 

Figure [12] shows the median values for the total IR lu- 
minosity (Ldust) as a function of flux in the PACS 160 and 
SPIRE 250 /xm bands at different redshifts. This plot illus- 
trates that, at low redshift, only a very few of the bright- 
est galaxies in our (relatively small volume) simulations 
would not even be considered LIRGs (Lir > 1O^^L0). Our 
simulations do not contain any ULIRGs (Lir > 1O^^L0) 
at low redshift. At ^ > 1, all galaxies above our nomi- 
nal GOODS- Herschel detection limit would be considered 
LIRGs or ULIRGs, at z ^ 2 all of these galaxies would 
qualify as very bright LIRGs or ULIRGs, and at 2; > 3 all 
of these galaxies are ULIRGs. 

At high redshift (2 < z < 4) our model predicts median 
SFRs that are a factor of two higher than those reported in 
"Lacey et al. (2010), as show in Fig. 11 As already noted, 
Lacey et al. (2010) assumed a top-heavy IMF in starbursts. 
With a top-heavy IMF, in a given star formation episode, 
a larger fraction of the mass goes into massive stars that 
can efficiently heat the dust. Thus, for a top-heavy IMF, 
the relationship between SFR and IR luminosity is shifted 
such that a given SFR results in higher IR luminosities. In 



the Lacey et al. (2010) model, the fraction of star formation 



occurring in the burst mode increases with redshift, so the 
average IMF with which stars are being formed shifts from 
being close to a solar neighbourhood IMF at the present 
day to being very top-heavy at high redshift. This helps to 
explain why our SFR vs 250 /xm flux relationship agrees with 



that of Lacey et al. (2010) in the local Universe (z < 0.25), 
but differs significantly at earlier cosmic epochs. In addition, 
the treatment of dust is also different in the [Lacey et al.| 
(2010) model. They used results from the RT+dust model 



calculations of Silva et al. ( 1998 ) 
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logio(6'i6o [mJy]) logio(5'250 [mJy]) 



Figure 8. Predicted median stellar mass vs. flux in the PACS 160 and SPIRE 250 fim bands for galaxies selected at different redshifts. 
The green dotted vertical line shows a representative flux limit of the GOODS-N observations. The error bars show the 16th and 84th 
percentiles. 




Iogio(^i60 [mJy]) logio(5'250 [mJy]) 

Figure 9. Predicted median dark matter halo mass (Mdm) as a function of IR flux in the PACS 160 and SPIRE 250 /im bands at different 
redshifts. The green dotted vertical line shows a representative flux limit for the GOODS-N observations. The error bars show the 16th 
and 84th percentiles. 



4.1 Sizes of High-redshift Late-type Galaxies 



For this sub-section we use our "high-redshift" SPIRE 250 
/xm selected sample of galaxies selected from our light- 
cones in the redshift range 2 < z < 4. We identify disk 
galaxies using the bulge to stellar mass ratio: galaxies with 
Mbuige/Mtotai < 0.4 are considered as late- type galaxies. We 
also removed all galaxies that have experienced one or more 
mergers within the last 500 Myr, because these galaxies have 
disturbed disks (and morphologies) and thus the size of the 
stellar disk is ill-defined. Figure [13] shows a Hess plot of the 
galaxy's disk size as a function of its stellar mass in 
physical units for this sample of high-redshift galaxies. For 
a comparison, we also plot the sizes and stellar masses from 
Cava et al. ( 2010J with green squares. The gray shading in 
both panels of Fig. [13] describes the number of simulated 



galaxies in a given two dimensional bin, while the red lines 
show the median and 16 and 84 percentiles. 

As in the real Universe, our model disks show a positive 
correlation between stellar mass and radial size. Very few 
galaxies are predicted to be as large in radial extent as the 



galaxies presented in Cava et aL] (2010) . However, the right 
panel of Figure [13] shows that when we select galaxies in the 
same redshift range with S'250 > 5 mJy, the predicted sizes 
are in good agreement with the observational results. The 
late-type galaxies in our 5*250 > 5 mJy high-redshift sample 
have average disk radii of ^ 2.2 kpc. This is significantly 
larger than the mean disk size (^ 0.9 kpc) of all late-type 
galaxies in the same redshift range (2 < z < 4). However, 
the small mean disk size in case of all late-types is driven 
by the large number of dwarf galaxies. But even if we use 
the same stellar mass cut that the SPIRE 250 /xm band 
selection induces, we find that the average size of the stellar 
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Figure 10. Predicted median cold gas mass vs flux in the PACS 160 and SPIRE 250 /im bands for galaxies selected at different redshifts. 
The green dotted vertical line shows a representative flux limit of the GOODS-N observations. The error bars show the 16th and 84th 
percentiles. 



3.0| ^ ■ ^ ^ ^ 1 3.0 




0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 

logio(5'i6o [mJy]) logio(5'250 [mJy]) 



Figure 11. Predicted median star formation rate as a function of IR flux in the PACS 160 and SPIRE 250 jim bands at different redshifts. 
The green dotted vertical line shows a representative flux limit for the GOODS-N observations. The error bars show the 16th and 84th 
percentiles. The hatched lines on the right hand side plot (LIO) are the predictions fromjLacey et alT](|2010|. 



disks is almost a factor of two smaller 1.3 kpc) than the 
mean disk size of the high-redshift sample. Consequently, 
in our models, the IR bright late-type galaxies that have 
not experienced a merger within the last 500 Myr are on 
average larger than galaxies that emit lower SPIRE 250 /xm 
fluxes. We speculate that this is because our models predict 
that larger galaxies are more gas rich, and we have already 
seen that there is a correlation between gas mass and IR 
luminosity. Our high-redshift sample of luminous IR galaxies 
also contains a few late-type galaxies whose stellar disks are 



> 7 kpc in size (for similar findings, see Rujopakarn et al. 



2011) 



4.2 Correlation with Merger Activity 

Ultra-luminous galaxies in the local Universe are widely be- 
lieved to be caused by major galaxy mergers (e.g. [Sanders 



Mirab el" 1996"; Coli na et al.|2001| [Farrah et al.|2003| [Dasyra 
et al.|2006b i Vaisanen et al.||2008 ). It has also been argued 
that ULIRGs will evolve into moderate-mass ellipticals (e.g. 



Genzel et al.|2001| [Tacconi et al. 12002) [Dasyra et al.|2006b| ) 



and/or AGN ( [Farrah et al. 2003 and references therein) 
after the merger or interaction. Merger driven star forma- 
tion has also been studied using hydrodynamic simulations 
e.g. [Jog k Solomon[[1992| [Mihos k Hernquist][1996[ [Cox 



et al. 2008| TVyssier, Chapon Bournaud|2010 ), which sup- 
port the picture that massive starbursts could be driven by 



mergers. It has been argued (e.g. Dasyra et al. 2006a) that 
the majority of ULIRGs are triggered by almost equal-mass 
major mergers and that less violent mergers of mass ratio 
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logio(5'i6o [mJy]) logio(^250 [mJy]) 

Figure 12. Predicted median total IR luminosity (Ldust) cis a function of IR flux in the PACS 160 and SPIRE 250 fim bands at different 
redshifts. The green dotted vertical line shows a representative flux limit for the GOODS-N observations. The error bars show the 16th 
and 84th percentiles. 
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Figure 13. Hess plot of disk sizes (i^disk) of the late-type galaxies 
as a function of stellar masses M^. All galaxies shown are in the 
redshift range 2 < z < A and have not experienced a merger 
within the last 500 Myr. The right-hand panel shows only late- 
type galaxies that are brighter than 5 mJy in the SPIRE 250 jim 
band. The red lines show the median and 16 and 84 percentiles for 
simulated galaxies. The observed late- type galaxies, shown with 
green squares, are from [Cava et al.| ( [20Tot . 



< 1 : 3 typically do not force enough gas into the nucleus to 
generate ULIRG luminosities. This is in agreement with sim- 
ulations (e.g. < 



& Bournaud 



Cox et al. ( 2008 ), but see also Teyssier, Chapon 



(2010) and references therein), although other 



factors (orbit, gas fraction, morphology) are also important. 
Even though some observations of IR luminous galaxies have 
failed to identify evidence for mergers, it has been argued 
that these highly obscured objects could be the end prod- 
ucts of collisions between gas rich galaxies in a final phase 
of merger evolution (e.g. [Auriere et aL| l996). Although this 
picture is well supported in the local Universe, the physical 
origin of ULIRGs at high redshift is less clear. For example. 



observations that suggest that high redshift galaxies may be 
able to achieve LIRG and even ULIRG luminosities with- 
out experiencing mergers. Our model assumes that mergers 
can trigger bursts of star formation, but also incorporates 
the high gas accretion rates at high redshift predicted by 
cosmological simulations. It is therefore quite interesting to 
investigate what our model predicts for the fraction of lumi- 
nous IR galaxies that are the result of mergers at high and 
low redshift. 

To study the fraction of IR bright galaxies that are pow- 
ered by mergers we record the time since the last merger 

^merger Or major merger Tmajormerger for each galaxy in 

our simulation. "Minor" mergers are defined as those with 
a mass raticQ greater than 1:10, and "major mergers" as 
those with mass ratio greater than 1:4. Our model predicts 
that about 84 (53) per cent of high-redshift galaxies with 
5*160 > 10 mJy have experienced a (major) merger dur- 
ing their lifetime. Note, however, that the fraction drops to 
about 34 percent if we concentrate on major mergers that 
have taken place less than 250 Myr ago, which are the merg- 
ers that are most likely to be causally linked to the star for- 
mation we are seeing. For the SPIRE 250 /xm band, the re- 
sults are very similar; about 86 (57) per cent of high-redshift 
galaxies with 5250 > 20 mJy have experienced a (major) 
merger. If we again concentrate only on recent mergers in 
which the major merging event took place less than 250 Myr 
ago, the merger fraction drops to ^ 42 per cent. Thus our 
model makes an interesting prediction: roughly every sec- 
ond high-redshift ULIRG can achieve ULIRG luminosities 
without a recent major merger. This prediction differs from 



the findings of Gonzalez et al. (2010), who used a semi- 
analytic galaxy formation model ("Baugh et al. 2005b) to 
study submillimetre galaxies (SMGs; ^sso > 5 mJy). Gon- 



Sturm et al. (2010) identified two galaxies in early Herschel 



[zalez et al.| ( |2010| finds that the majority of SMGs seen in 
current sub-mm surveys are starbursts triggered by major 



^ Here, the masses used in the ratio include stars, gas, and the 
dark matter within the optical radius of the galaxy, as in S08 
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Figure 14. Merger fraction as a function of PACS 160 fim flux 
at different redshifts. The upper-left plot shows the recent major 
mergers, while the lower-left shows recent minor mergers (see text 
for deflnitions). The upper right panel shows the merger fractions 
as a function of PACS 160 /im flux for all mergers within the 
past 500 Myr, and the lower right panel shows the fraction of 
galaxies that have never merged. The vertical dotted line shows 
a representative flux limit for GOOT)S- Herschel (5 mJy). 
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Figure 15. The fraction of galaxies that have experienced a 
merger and that have never merged with another galaxy as a 
function of far- to mid-IR flux ratio. The flux ratio is the ratio 
between the PACS 160 and Spitzer IRAC 4.5 /xm bands. The 
merger events have been divided into recent or young (< 250 
Myr) and less recent or old (250 < Tmerger < 500 Myr) and to 
minor and major incidents and are colour coded as given by the 
legend. 



and minor mergers involving gas-rich disk galaxies (see also 
Baugh et al.||2005b| ). 

We present these results in more detail in Figure |14| 
which shows the evolution of merger fractions as a function 
of PACS 160 /im flux at different redshifts. The upper-left 
plot shows that the merger fraction of recent (the latest 
merger took place less than or equal to 250 Myr ago) ma- 
jor mergers grows as a function of PACS 160 /xm flux. This 
is true independent of studied redshift, although the steep- 
est rise is noted at z = 1 and 4. The plot shows that the 
likelihood of the most luminous ULIRGs being powered by 
a merger is higher than for those that are less luminous. 
A similar trend can also be seen for minor mergers. This 
shows that under the right conditions, minor mergers may 



also trigger (U)LIRG-like activity. We note that Fig. 14 and 
the related results are qualitatively similar if instead of the 
PACS 160 /xm flux we use the SPIRE 250 /xm band. 

We can use our models to attempt to identify obser- 
vational quantities that could provide interesting indicators 
of starburst and merger activity for observations that are 
available or will become available soon. We wish to exploit 
the availability of extensive mult i- wavelength data, from the 
rest UV to the FIR with Herschel, in fields such as GOODS. 
We focus on the redshift range 2 < z < 4 for these experi- 
ments, and use the high redshift sample from our lightcones. 
It is reasonable to expect that starbursts will have high val- 
ues of specific star formation rate (SSFR), M^^ . In a 
relatively narrow redshift range, we can use the ratio of a 
FIR fiux to a Spitzer IRAC band as a proxy for SSFR. 

Figure [15] shows the fraction of different types of merg- 
ers and the fraction of galaxies which have never merged as 
a function of far- to mid-IR fiux ratio (colour). We choose 
to use the PACS 160 /xm band as our far-IR band, because 
it is relatively close to the peak of the IR fiux distribution 
in the studied redshift range (2 < 2; < 4), but at the same 



time it suffers less from confusion than, e.g., the SPIRE 250 
/im band. The mid-IR band chosen is the channel two of the 
Spitzer IRAC at 4.5 /xm, which should provide a proxy for 
the stellar mass of a galaxy. Fig. [15] shows that our models 
predict that the probability that a galaxy has experienced a 
recent merger is a strong function of this far- to mid-IR fiux 
ratio. For the most extreme values (Siqo S^l ^ 1300) more 
than half of the galaxies have experienced a major merger 
within the last 250 Myr. We note that the results remain 
if we instead use other bands such as 5'250 Sg'^ as the fiux 
ratio. It will be interesting to test whether these IR colors 
are correlated with signatures of recent merger activity in 
observed galaxies. 

We can also explore whether the use of a two colour di- 
agram can provide further information. We expect that the 
rest-frame UV colours might correlate with the age of the 
starburst, and thus with the time since the merger event. 
We consider the filters of the Hubble Space Telescope's Ad- 
vanced Camera for Surveys, because such data are available 
for GOODS. In Figure [16] we combine the information on 
UV-optical colour and far- to mid-IR fiux ratio in a colour- 
colour plot. Again, we choose to use the PACS 160 /xm band 
and IRAC's channel 2 (4.5 /xm), because these channels allow 
the most robust comparison with observations of galaxies in 
the GOODS field. However, the results would be qualita- 
tively similar with any IRAC channel or with the SPIRE 
250 /xm band. Fig. [16] shows that in our models, the most 
recent mergers are all located in the lower right region of 
each subplot, that is the region with the most extreme far- 
to mid-IR fiux ratio and the bluest UV-optical colour. 

The black contours indicate where the galaxies which 
have never merged are located and have been drawn using 
two dimensional Gaussian kernel density estimator. Note, 
however, that this region is also populated by recent mergers 
so we cannot conclude that the observed galaxies in this 
region have not experienced a recent merger. However, it is 
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interesting that there is a region of the diagram (shown by 
a green dotted hne) in which all galaxies have experienced 
a recent merger. The functional form of the green line is as 
follows: 

{F775W - FSmp) AB = A log^o ^ + B , (3) 

04.5 

where A = 1 and B — —2.4. Note, however, that this simple 
selection does not seem to be able to separate recent major 
mergers from less extreme mergers. If instead of the PACS 
we had used the SPIRE 250 /im band, the above equation 
would remain the same, but with B = —2.7. 

It will be very interesting to test these predictions using 
surveys such as GOODS- Herschel, where multi-wavelength 
photometry as well as high resolution Hubble imaging are 
available. We plan to pursue this in a future work (Niemi et 
al. in prep). 



5 SUMMARY AND CONCLUSIONS 

In this work we made use of the semi-analytic models devel- 
oped by Sll, in which attenuation and re-emission of light 
by dust is modelled with a simple analytic approach, to 
make predictions for the physical and observable properties 
of galaxies that may be detected in deep fields by the Her- 
schel Space Observatory. We used the semi-analytic models 
to generate mock lightcones resembling the GOODS survey 
but with 100 times the volume. 

We first examined the differential number counts for 
galaxies in the PACS and SPIRE bands. We found very good 
agreement with the observed number counts derived from 
Herschel PACS data, for the overall counts and for galaxies 
binned by redshift at z < 2. At 2; > 2 our model underpre- 
dicts the number of more luminous galaxies in the PACS 
100 and 160 /xm bands. The agreement of our model predic- 
tions with the observed counts deteriorates as one moves to 
longer wavelengths, with increasingly poor agreement seen 
as one moves from PACS to SPIRE 250 fim and then to the 
longer SPIRE bands. At 250 /xm and longer wavelengths, the 
models tend to underpredict the number counts of interme- 
diate flux galaxies (10 ^ S250 ^ 100 mJy), and overpredict 
the counts of brighter galaxies (5*250 ^ 100 mJy). We pre- 
sented predictions using three different sets of empirical dust 
emission templates, and found interesting differences. How- 
ever, it does not appear that any single simple modification 
to the dust templates could solve the discrepancies. 

We discussed possible reasons for these discrepancies. 
These include uncertainties due to cosmic variance, possible 
evolution or scatter in the dust emission templates, a possi- 
ble contribution to the IR light from obscured AGN, and the 
possibility that many of the bright sources are lensed. It is 
also possible that some of the physical processes in our model 
need revision, or that the stellar IMF was different in early 



galaxies, as has been suggested by Baugh et al. (2005b) 



However, our models produce reasonable agreement with the 
luminosity functions at rest-UV through NIR wavelengths, 
and with the observationally derived stellar mass functions 



and SFR functions at high redshift (see Sll; Fontanot et al. 
2009a and references therein). Moreover, the top-heavy IMF 
model of 'Baugh et al.' {'2005b ) that was developed to ex- 
plain the observed sub-mm population at 850 /xm does not 
appear to reproduce the SPIRE number counts particularly 



well (Cle ments et al.|2010| Lacey et al.|2010 ), and does not 
produce enough galaxies with large stellar masses at high 
redshift ( Swinbank et al.|2008 ). We also discuss the impact 
of flux boosting and blending on the observationally derived 
Herschel counts. This would explain why the discrepancy 
becomes larger with increasing wavelength, and might also 
be a larger effect for high redshift galaxies, which are com- 
pact and strongly clustered. A quick visual inspection of 
HST images for high redshift galaxies makes it clear that 
many, if not most, Herschel sources probably have contri- 
butions from multiple galaxies. Moreover, according to the 
simulations carried out by Rigby et al. (2011), more than 
half (^ 57 per cent) of sources detected at > 5cr at 500 /im 
show a flux boosting by a factor > 1.5 and more than ev- 
ery fourth (^ 27 per cent) by a factor > 2. This is clearly 
a very serious issue when comparing theoretical predictions 
with these observations. We intend to use our simulations to 
create mock Herschel images, carry out the source detection 
procedure on these, and study the predicted effects of blend- 
ing on the measured fluxes and galaxy counts. In addition, 
the Herschel team is working hard to better understand and 
overcome this problem. 

With the resulting caveat that one probably cannot in- 
terpret our predicted Herschel fluxes as corresponding in a 
one-to-one fashion with directly observed Herschel fluxes (as 
currently available), but rather as intrinsic fluxes in the ab- 
sence of blending, we used our model to make predictions for 
physical properties of Herschel detected galaxies. We showed 
the predicted relationships between PACS 160 and SPIRE 
250 /xm fluxes and stellar mass, halo mass, cold gas mass, 
SFR, and total IR luminosity as various redshifts. Our mod- 
els show fairly strong trends between both ^leo and 5250 at 
a given redshift and all of these quantities except halo mass, 
which is quite flat as a function of both 5i6o and ^250 at 
higher redshifts. This may have interesting implications for 
the clustering properties of these sources (for early observa- 
tional results in GOODS South, see e.g. Magli occhetti et al."] 



2011). Our models predict that the galaxies that are likely 



to be detected in deep Herschel surveys such as GOODS- 
Herschel {Sieo or ^250 > 5 mJy) at high redshift {z > 2) 
have fairly large stellar masses (M^ 6 x 10^° -3 x lO^^M©) 
and reside in fairly massive dark matter halos ( ^ 10^^ M©). 
The latter is in good agreement with the minimum mass de- 
rived by Amblard et al. (2011) based on the observed clus- 



tering of SPIRE sources. 

The simulated IR luminous galaxies were also found to 
be gas rich with cold gas masses ranging from ~ 10^ up 
to - 4 X lO^^M©. The bulk of the high redshift IR bright 
galaxies were found to contain ^ 6 x lO^^M© of cold gas. 
These predictions will soon be able to be tested with radio 
and mm observations such as those that will be obtainable 
with ALMA. 

Our predicted stellar, cold gas, and halo masses at a 
given 250 fim flux are significantly higher than the corre- 
sponding predictions from the model of Lacey et al. (2010), 



which assumes a top-heavy IMF in starbursts. We show a 
quantitative comparison between our predicted SFR vs. ^250 
relation and that of Lacey et al. (2010), finding that al- 



though their predictions are similar to ours for low redshift 
galaxies, at higher redshift {z ^ 0.5), their predicted SFR 
are lower than ours (at a given FIR flux) by about a factor 
of two. 
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Simulated Galaxies; Sim > 5 mJy 




Figure 16. Rest-frame UV-optical colour as a function of far- to mid-IR flux. The flux ratio is the ratio between the PACS 160 
and Spitzer IRAC 4.5 /im bands. The model galaxies have been colour coded by the time since the last merger in Myr. The green 
dotted line marks the region in which all simulated galaxies are recent mergers. 



We investigated the radial sizes of high redshift (2 < 
z < 4) late-type galaxies selected in the SPIRE 250 /xm 
band and compared our predictions with the observations 



recently presented by |Cava et al. (2010). We found that the 
SPIRE detectable galaxies are on average significantly larger 
than IR faint galaxies with comparable stellar masses. The 
mean disk size of IR luminous late-type galaxies was found 
to be in good agreement with the observations, although 
the statistics are extremely limited. We speculate that this 
is because larger galaxies are more gas rich and therefore 
more likely to be IR bright. 

Finally, we investigated our model predictions for the 
importance of merger-driven starbursts in producing IR lu- 
minous galaxies at different redshift s. Although there is 



clear evidence for a connection between merger activity and 
LIRGs and ULIRGs in the local Universe, the situation at 
high redshift is less clear. However, it has been suggested 
that an increasing fraction of LIRGs and ULIRGs at high 
redshift might not have a merger origin. We present quanti- 
tative predictions for the fraction of galaxies that have ex- 
perienced a major or minor merger within a given timescale 
as a function of redshift and 160 fim flux. We find a fairly 
strong trend between the 160 fim flux and the probabil- 
ity that a galaxy has had a recent merger, indicating that 
brighter galaxies are more likely to be merger driven. How- 
ever, we find the interesting result that in our models, a sig- 
nificant fraction (half or more) of IR-luminous galaxies at 
high redshift {z > 2) have not experienced a recent merger. 
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This implies that the high gas accretion rates and efficient 
feeding via cold flows predicted by cosmological simulations 
at high redshift can fuel a significant fraction of the galaxies 
detected by Herschel. This appears consistent with prelimi- 
nary observational results (e.g. [Sturm et al.||2010| , but will 
be interesting to study and quantify in more detail in the 
near future. 

In order to pursue this question further, we used our 
model to try to identify some observational photometric sig- 
natures associated with galaxies that have experienced re- 
cent mergers. We showed that there is a strong correlation 
between the mid-to-IR colour Siqq/Sa.^ (PACS 160 /xm to 
Spitzer IRAC channel 2 ratio) and the probability that a 
galaxy has had a recent merger. We also investigated a two 
color diagnostic using rest-frame UV colour and the 160 /xm 
to IRAC 4.5 /im ratio. We found that recent mergers occupy 
a characteristic region of this diagram, suggesting another 
interesting observational tool to identify probable recent 
mergers. In a work in progress, we are making use of HST 
imaging in the GOODS- Herschel survey to study the mor- 
phologies of galaxies selected using these multi-wavelength 
diagnostics, in order to test our model predictions and the 
importance of merger-induced star formation at high red- 
shift (Niemi et al., in prep). 
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